Mott Insulator - Superfluid Transitions in a Two Band Model at Finite Temperature 

and Possible Application to Supersolid "^He 
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We study Mott insulator - superfluid transition in a two-band boson Hubbard model, which can 
be mapped onto a spin-1/2 XY model with spins coupled to an additional Ising degree of freedom. 
By using a modified mean field theory that include the effects of phase fluctuations, we show that 
the transition is first order at both zero and finite temperatures. On the Mott insulator side, there 
may be reentrance in phase transition. These features are consequences of the underlying transition 
between competing defect poor and defect rich phases. The relevance of the model and our results 
to supersolid "^He and cold bosonic atoms in optical lattices are discussed. 



I. INTRODUCTION 

The Mott insulator- Superfluid (MFSF) transition is 
one of the striking phenomena arising from the many 
body physics of correlated boson systems on a lattice^. 
The physics of this transition is believed to be relevant 
to not only systems with bosons as "elementary parti- 
cles" , but also condensed matter systems whose low en- 
ergy physics is governed by bosonic degrees of freedom, 
for example, an array of quantum Josephson junctions. 
Recently, there has been a renewed interest in the MF 
SF transition due to its relevance to two experimental 
systems. The first is supersolidity of solid ''He^i^ if one 
views the normal solid to supersolid transition as a Bose 
condensation transition that occurs due to delocalization 
of the ''He atoms without melting the underlying crys- 
talline lattice. The second is ultracold bosonic atoms in 
an optical lattice^. 

The commonly accepted paradigm for the MFSF tran- 
sition is based on the one-band boson Hubbard model. 
In that model, the MI-SF transition is a continuous one 
with the Mott-Hubbard gap and the condensate density 
increasing continuously from on the MI and the SF sides 
of the transition respectively. Recently, Dai, Ma, and 
Zhang (DMZ)^ proposed a two-band Hubbard hamilto- 
nian as a model for solid "^He, and showed using a single- 
site mean field theory that the Ml-supersolid transition 
can be first order. In this paper, we seek to investigate 
the validity of this claim beyond mean field theory by 
including the effects of phase fluctuations of the conden- 
sate. Specifically, we consider the following spin S — 
XY model on a lattice, where the spins are coupled to an 
additional Ising degree of freedom 
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Here the Ising variable fii — 0,1 and and S~ are the 
raising and lowering operators. As a model of magnetism, 
this Hamiltonan would describe a XY spin system in a 
transverse magnetic field h, with = 0, 1 representing 
the absence (presence) of a spin on site i. In this case A 



plays the role of a chemical potential for the spin. For 
example, H can act as a model for an electronic insulator 
with two on-site electronic configurations, one magnetic 
and one non-magnetic separated by energy A, or it can 
describe a binary alloy of magnetic and non-magnetic 
ions. However, for this work, we are more interested in 
H as SL model of boson insulator-superfluid transition. 
Below we will discuss the relevance of this model to the 
boson two-band Hubbard model and to supersolid "^He 
and trapped bosons in an optical lattice. 



A. Supersolid He 

The possibility of supersolidity in ^Hc, which refers to 
the coexistence of crystalline ordering and superfluidity, 
was proposed some years ago by Andreev and Lifshitz », 
Chester— , and Leggeti*^. The interest in such an unique 
state of matter is recently rekindled following the report 
by Kim and Chan2»^ of the observation of non-classical 
rotational inertia (NCRI) in solid ''He confined inside 
porous media and in bulk solid ''He at low temperature. 
NCRI has since been confirmed by other groups^ i^^'^"' , 
although it is still debatable as to whether this is a 
bulk equilibrium effect, or a result of non-equilibrium 
defects^ii^. Thus far, experiments have failed to detect 
direct superflow or " superconductivity" ^'^i^^ . Therefore, 
it remains controversial as to whether ^He is a supersolid 
at low temperature. Nevertheless, these recent experi- 
ments provide the impetus for better theoretical under- 
standing of the mechanism for bulk equilibrium super- 
solidity. In their seminal paper, Andreev and LifshitZi^ 
proposed that since ^He is a quantum solid with large 
zero point motion, a finite density of defects may be 
present even in the ground state. Amongst these de- 
fects, the most promising candidate is zero point vacan- 
cies, whose motion they argued will be wave-like at low 
temperature and so will Bose condense. The Bose con- 
densation of the vacancies can then lead to superfluidity 
without destroying the underlying crystalline ordering. 
From the point of view of the presence of off-diagonal 
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long ranged order (ODLRO) in Jastrow wavefunctions 
describing solids, Chester- also conjectured zero point va- 
cancies as the mechanism for supersolidity. If vacancies 
are present, the solid will be incommensurate, with the 
number of He atoms different from the number of lattice 
sites. Recently, Anderson, Brinkman and Huse (ABH)i^ 
pointed to the correction to the specific heat as evi- 
dence of incommensurability in the ground state. 

While appealing, the idea of zero point vacancies is 
subject to stringent constraint from both experiments 
and several computational calculations, showing vacan- 
cies as activated from the commensurate solid. The ac- 
tivation energy is found to be about 10 A' from X-ray 
diffraction^*' and about 15 K in simulations^^. This is at 
odds with the Andreev and Lishitz's proposal that the 
commensurate ground state is unstable with respect to 
spontaneous generation of vacancies. Recently, Dai, Ma, 
and Zhang (DMZ)Si proposed a possible solution to the 
quandary by including a lower energy defect that they 
called the exciton, which is an on-site bound state of a 
vacancy and an interstitial. Physically, this defect corre- 
sponds to a broadening distortion of the local wavefunc- 
tion. For example, if the defect free solid state is given by 
the (unsymmeterized) Hartree-Nosanow wavefunctioni- 



where (f>^ is a localized wavefunction, then the state with 
one exciton defect on the site j can be written as 

where is also localized, but less so than 0^. The key 
idea is that an atom in the state 0;, can tunnel more ef- 
fectively into a neighboring vacancy site (and vice versa) 
than one in the state (j)^ due to the wavefunction broad- 
ening. As a result, even though vacancies are activated 
in a defect free background, they can be spontaneous 
generated in an exciton background. At low tempera- 
ture, they will then Bose condense and the condensa- 
tion energy may overcome the energy cost of creating 
exciton defects. In their theory, the commensurate solid 
would then be a metastable state, while the true ground 
state will be an incommensurate one with finite densi- 
ties of vacancies and excitons. The metastability of the 
commensurate solid has also been proposed by Andeson, 
Brinkman, and Huse^^. 

As a model for their theory, DMZ consider a two-band 
Hubbard model, in which the defect free (DF) state is 
considered the "vacuum", and the boson operators are 
operators that create and destroy defects. In this paper, 
it is more convenient to take the vacuum as the physical 
vacuum for "^He atoms, in which case the DMZ Hamilto- 



nian becomes 

H = (egg+ai + ebhjh + Uhgifn,^ (1) 

i 

-^{iaatcLj + Ubfbj + h.C.) 
(il) 

Here, the strong repulsion between '^He atoms in close 
proximity is modeled by taking a and &, field operators 
for "'He atom in the a and b states, as hard-core boson 
operators. U is the repulsion between a ^He atom occu- 
pying (j)g and one occupying (/(^ on the same site, and is of 
the order of the local interstitial activation energy. Of the 
various energies in the problem, this is by far the largest 
50 K) , and for the purpose of this work we will con- 
sider f/ — > oo, so each site can at most be singly occupied 
by a *He atom. In other words, we neglect the possibility 
of interstitials. ta and tf, denote the hopping amplitudes 
of the tunnelling process between neighboring sites from 
a to a-state and from b to &-state respectively {{ij) indi- 
cates nearest neighbors). The 5-boson has higher on-site 
energy (sf, > e^), but also a higher hopping amplitude 
{tb > ta). The hard core nature of the bosons implies 
that hopping is only possible between an occupied site 
and a neighboring vacancy. 

In the torsional oscillator experiment setup used to 
observe NCRI, the density of ^He atoms is held fixed. 
However, because the lattice constant of the solid is de- 
termined by minimizing the free energy, the number of 
bosons per site {fia) + {fib) is not externally imposed. 
Within this model, if {fia) + {fib) = 1, the '*He atoms 
form a commensurate solid. If further {fib) = 0, the com- 
mensurate solid is defect free (DF), while if {ha) — 0, it 
is an exciton solid. However, if {ha) + {hb) < 1, then 
the solid is an incommensurate one with a finite vacancy 
density. The excitation energy of a single vacancy from 
the DF state is |£a| — zta, where z is the coordination 
number. The Andrccv-Lifshitz vacancy mechanism for 
supersolidity would require this to be negative, which 
evidently is not the case. On the other hand, the insta- 
bility criteria for generating a vacancy above the exciton 
solid is \ea\ — ztb, which is easier to satisfy. Vacancies 
can then Bose condense above the exciton background 
and the condensation energy gain may be sufficient to 
overcome the exciton energy A — Sb — Sa- Should this be 
the case, the DF state will be metastable. In their theory, 
DMZ identify the DF state as the normal solid at T = 0, 
and the exciton state with Bose condensation of vacancies 
as the supersolid. Using a single-site mean field theory, 
DMZ showed that as parameters in the mode are tuned 
(corresponding experimentally to for example changing 
the pressure), there is a transition at T = from the 
DF normal solid to the supersolid state described above. 
Since the DF state is metastable, this transition is natu- 
rally first order, and involves a jump not only in the Bose 
condensed amplitude, but also in the densities of both va- 
cancies and exciton defects. Hence the supersolid transi- 
tion is accompanied by a commensurate-incommensurate 
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transition, as well as a change in local density profile of 
the ^He atoms. Because the normal solid is commen- 
surate, it explains why experimental measurements per- 
formed on the normal solid do not show an appreciable 
density of vacancies. 

Because the results of DMZ are based on a single-site 
mean field theory, it is of interest to examine if they hold 
up against the effects of both quantum and thermal fluc- 
tuations. At low temperature, the dominating effects 
should be due to phase fluctuations of the condensate, 
since they correspond to gapless modes. Experimentally, 
at the lowest temperature, NCRI in solid ^He is observed 
even at the highest pressure achievable in the laboratory. 
Thus, the T = first order supersolid-normal solid tran- 
sition predicted by DMZ cannot be tested at present. 
It is therefore important to determine the nature of the 
supersolid-normal solid transition at finite temperature 
within DMZ's theory to compare to experiments. We 
confirm that when phase fiuctuations are included, the 
finite temperature transition remains first order. More- 
over, we find that fluctuations give rise to a 'reentrance' 
phenomenon. 



B. Trapped Bosons in Optical Lattice 

For solid ^He, the lattice constant is self-adjusted to 
minimize the free energy. We can also consider a system 
of bosons in the presence of an external periodic poten- 
tial, in which case the lattice constant will be externally 
imposed. An interesting realization of such a system has 
been achieved recently in ultracold trapped bosons in an 
optical lattice^ '^^i^° . By superimposing counter propa- 
gating laser beams of wavelength A in different direc- 
tions, an effective periodic potential in one, two, or three 
dimensions can be produced. In 3-D, the potential has 
the form 

V{x, y, z) — Vox cos^lkx) + Voy cos^{ky) + Vqz cos^(fcz) 

where k = 2t:/X is the wavenumber, and the strengths 
of the potential in the three directions are proportional 
to the laser intensity in each direction and can ther- 
fore be tuned. For a single particle, this periodic po- 
tential gives rise to the usual energy bands in the flrst 
Brillouin zone. For Vq large compared to the so-called 
recoil energy En = fi^fc^/2m, the lower bands can be 
viewed as tight-binding bands arising from a periodic 
array of harmonic wells. For deep potential of cubic 
symmetry {Vqx = Voy = Vqz = Vq), the Wannier basis 
for these bands can to a good approximation be con- 
structed from the cigenstates of the spherically symmet- 
ric harmonic oscillator potential Vhar{r) — Vok^r^. The 
bands can then be denoted by the quantum numbers 
{nx,ny,nz) of these eignestates, and their eigenenergies 
(rix + riy + riz + 1/2) Hloq, where Hujq = \/AVqEj^, pro- 
vide approximate values for the band gaps. It is conve- 
nient to consider the optical lattice experiments as being 



performed at fixed chemical potential. In the case where 
HujQ is much larger than other relevant energy scales, all 
the particles will occupy the lowest band, and the sys- 
tem can be modeled by a single-band boson Hubbard 
mode&: 

i i 
{ij) i 

with fi as the chemical potential. This is a well-studied 
model, and it has been established that with increasing 
ta/U, there is a quantum phase transition from a Mott 
insulator state with integer filling per site to a superfluid 
(SF) state. The SF state may have commensurate or in- 
commensurate filling depending on fj,. Trapped bosons in 
the optical lattice provides an elegant realization of this 
model which allows the transition to be studied system- 
atically in experiments. 

Alternatively, we can consider the limit where the Hub- 
bard U is the dominating energy. This can be achieved 
experimentally by tuning close to the Feshbach reso- 
nance. In this case the Mott-Hubbard gap prevents any 
double or higher occupancy Ui < 1, and the bosons are 
hard-core. If we restrict to the lowest band still, then 
the ground state is an insulator for /i < Eq — zta (vac- 
uum) and ^j, > Sa + zta (one boson on each site), while 
if /i lies between these two limits, the system is a su- 
perfluid with (rii) changing continuously from to 1 as 
fi increases. Unlike the MI-SF transition with changing 
ta/U, these insulator-SF transitions caused by changing 
fi in the hard-core model are just by-products of density 
transitions. 

For ^ > Ea + zta, one has essentially a filled band 
insulator if we restrict to the lowest band. However, if the 
higher bands are considered, then it is possible for bosons 
to delocalize and form a SF. We emphasize from the start 
that this is not simply due to partially occupying the 
higher bands while keeping the lowest band filled or even 
transferring some bosons from the lowest band to a higher 
band, as the hard core condition would still forbid any 
boson motion if every site is singly occupied independent 
of which band it is in. Instead the boson delocalization 
can only occur by introducing vacancies. 

We consider the case where in addition to the low- 
est band, we include one higher band, with both strong 
enough intraband and interband on-site repulsion be- 
tween bosons to make them hard-core. One thus arrives 
at a two-band boson model 

H ^ Y (eaa^ai + Ebbth^ - ^^itaajaj + tbbjh.j + h.c.) 

i 

with the constraint Uai + nu < 1. This Hamiltonian is of 
the same form as the DMZ 2-band model for superolid 
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except for the presence of the chemical potential term, 
which can simply be absorbed into a redefinition of 
and El, . The d and b operators are for the lowest and 
the higher bands respectively. In analogy to the super- 
solid model, we will refer to a boson in the 6-band as an 
exciton. Because the Wannier state for the higher band 
is less localized, the hopping amplitude tb will have a 
bigger magnitude than ta if the lattice constant A/2 is 
sufficiently large. The ratio of the hopping amplitudes 
can be estimated by looking at the ratio of overlaps be- 
tween states on neighboring sites. Using this estimate, we 
calculate the quantity tmn/too^ where tmn is the proba- 
bility density of the boson tunneling from eigenstate n 
of a well into eigenstate to of a neighboring well. In Fig. 
[IJa), we plot too a function of the Gaussian half- width 
of the harmonic oscillator. The dependence of the ratio 
tmm/too on the distance between wells, A/2, is shown in 
Fig. Wljo) for TO = 0, 1, 2 in one dimension. If we identify 
a as the n = state, and b as either the n = 1 or n = 2 
state, then we can see clearly that tj, can become sig- 
nificantly larger than ta as X/2xq increases, where xq is 
the Gaussian half- width of the n = harmonic oscillator 
eigenstate. 

Instead of the lowest and a higher band of a single type 
of bosons, the two-band model also acts as the model for 
the case of two different types of bosons (a, b) each re- 
stricted to its respective lowest band in the optical lat- 
tice. For example, we may have bosons of two different 
masses (to^ > mi,). If we assume the laser beams pro- 
duce the same periodic potential on them, then we have 

A = Eb - Ea = fkoob- fiwoa, whcrc tujQa.b cc ^J^^ are 
the harmonic oscillator frequencies for mass to^ and to^ 
respectively. We also have tb > ta, due to larger zero 
point motion of the &— boson. However, because the 
tunneling amplitude ^qq^ is exponentially sensitive to the 
gaussian half-width of the harmonic oscillator wavefunc- 
tion, we expect the difference in tunneling amplitudes to 
be the more significant effect. 

Before we discuss our calculation for the two-band 
model, we should clarify the terminology used below. 
In both solid ^He and trapped bosons in optical lattice, 
Galilean invariance is violated. In the optical lattice case, 
the lack of translational invariance is externally imposed, 
so when bosons Bose condense, we consider it to be a su- 
perfluid, not a supersolid. In the case of solid "^He, the 
translational invariance is spontaneously broken by the 
solid, and if Bose condensation occurs without the solid 
melting, we have a supersolid, to be distinguished from 
the superfluid, which refers to the situation in the liquid. 
However, since we are not interested in the melting tran- 
sition, and since the two-band model is applied to both 
solid "'He and optical lattice, we will simply refer to the 
Bose condensed phase as the SF phase in both systems 
for convenience in what follows. 




W2x„ W2x„ 



FIG. 1: (a) The overlap too of the lowest energy states and (b) 
tmm {m = 0, 1,2) scaled by too, as functions of the distance 
A/2 between two neighboring centers in the unit of Gaus- 
sian half-width xo for a ID harmonic oscillator in estimat- 
ing the tunnelling amplitudes of bosons. We use tmn (a) = 
I!!°aa V'm (^) fn ~ , where (a; — a) denotes the eigen- 
state of the n*'' energy level of the harmonic oscillator with 
center x = a. 



II. MODIFIED MEAN FIELD THEORY 

Our primary goal is to analyze the above 2-band bo- 
son Hubbard model to include phase fluctuations. Since 
we do not expect the precise details of the individual pa- 
rameters in the Hamiltonian other than the gross features 
described above to be crucial, we take ia = for conve- 
nience. The boson model can then be mapped onto a 
spin S" = 1/2 XY model coupled to an Ising degree of 
freedom as follows. 

H = ^(ea-M) + ^(^A+^U,-/i^^,,n,(2) 

i i ^ i 

(v) 

with h = fi — {ea + A). The mapping between the original 
boson operators and the new set of operators n and S is 

|0) = |0) , nSt |0) = b+ |0) , h,a+ |0) = (3) 

S■=btb^-l, St = bt, S-=h (4) 

where |0) is the empty site, i.e. vacancy state. The Ising 
operator hi has eigenvalues and 1. For bosons, it acts 
as the defect (vacancy or exciton) occupation number 
operator. The spin 1/2 operator Si acts on the Hilbert 
space of the |0) (vacancy) and bf |0) (exciton) doublet, 
which is assumed to be closer to each other in energy 
than to that of |0) . 
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Because we take ta = 0, the a-bosons cannot Bose 
condense, and the Bose condensation order parameter, 

taken to be real, is given by (hi^ — (b^'^ — (^Six^ . In 
this representation, the single site MFT of DMZ would 
correspond to approximating S~ hifij by 

+n. (5+) (57) (n,) + h, (S+) (57) (n,) 

If we apply the single-site mean field theory to this model, 
we simply recover the results of DMZ for the choice of pa- 
rameters used here. Such MFT of course neglects fluctua- 
tions completely, and it is legitimate to question if the re- 
sults remain valid when fluctuations are included. At low 
temperature, the dominating fluctuations in three dimen- 
sions should be those from gapless excitations, which are 
the phase modes of the Bose condensate, or spin waves 
in the spin language. In order to include the effects of 
phase fluctusations, we modify the mean field theory as 
follows 

SiSjhihj = l^SiSj'^ {{hi) hj + hi (hj)) (5) 
+SiSj {hj) - 2 l^S^Sj'^ {hi) {hj) 

That is, correlations between the Ising degrees of free- 
dom on different sites and between the Ising and spin 
degrees of freedom on same or different sites are still ig- 
nored. However, correlations between the spin degrees 
of freedom on different sites will be included to allow for 
spin wave physics. The modified mean field Hamiltonian 
can now be written as 



are parameters independent of site number and to be 
determined self-consistently. It can be seen that the 
spin part Hamiltonian Hs is a ferromagnetic XY model 
in a transverse field with exchange coupling J . 

Hn is a single-site Hamiltonian and {hi) is easily cal- 
culated. The self-consistent equation for n is then given 
by 

where 

Ei = A + ^-hM - 2tbZhB 

The factor i in the denominator arises from the Ui = 
state being a singlet state while the rii = 1 state is a 
doublet. 

The most important feature of Hs is its continuous XY 
symmetry which gives rise to gapless excitations (Gold- 
stone modes) in the ordered state. This physics will not 
be affected by the transverse field h provided is not too 
big. Thus, as further simplification, we consider in this 
work the zero field case, or in other words the case where 
the 5-state and the vacancy are degenerate. Setting 

= 0, the three parts of the MF Hamiltonian becomes 

Hs = -jY.{s^iS^ + Sfsf) (10) 
H„ = {A - 2UZhB)Y,fl^ (11) 



C [{Ea - m) + 'iUZn^B] (12) 



where 



Hs 



-hhY,st-JY.(s^s^ + sy,sy 

i (ij) 



Hn= f A + - - hM - 2tbZhB j^n^ (7) 



I 






and Z is 


M 


- (^^^> 


B 




n 


= {h^) . 



(6) 



C = ^{(£a- At) +/inM + 24^^25} (8) 



A. The Modified Spin- Wave Method 

After the modified mean field decoupling, Hg remains 
a many body Hamiltonian. Since our main goal is to 
include the effects of the condensate phase fluctuations, 
which in the spin language implies fluctuations of spin 
waves, it is natural to use the spin wave approximation. 
The spin wave theory is known to be reliable in 3D, where 
long range order persists at low temperature. However, 
it is seen that the coupling J is proportional to n^, which 
is itself temperature dependent, and as we will see, will 
imply that we need to address Hs both below and above 
the spin ordering temperature. One possible method that 
can be used is the Schwinger boson MFT, but it has been 
shown that this method is problematic for temperature 
comparable or larger than J.^^ Instead, we use a modified 
spin wave method similar to the approaches introduced 
by Takahasi22i2i and by Tang and Hirsch— . 
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The calculation will be performed on the 3D cubic lat- 
tice. For the case of bosons on the optical lattice, this is 
the actual lattice that has been studied experimentally. 
For solid ''He, the lattice is hep, but since Hg is ferromag- 
netic, the physics does not depend on the precise lattice 
structure beyond the co-ordination number Z. To per- 
form the spin wave theory, it is convenient^^ to globally 
rotate the spins about the spin y-axis and rewrite Hg as 

IS 

where the P's are spin operators in the rotated frame. 
The classical ground state can then be taken as having 
all the spins pointing in the +z direction. Next, the 
Holstein-Primakoff transformation is defined 

P+ = ^2S~ctciCi 

Pr ^ 4^2S-cici (13) 

Pt = S~c+ci 

The linearized spin wave Hamiltonian Hsw is obtained 
by expanding the square roots and keeping to quadratic 
order in the bosonic c operators. This quadratic Hamilto- 
nian can then be solved using Bogiiubov transformation 

dk = cosh^fcCfe-sinh^fcC^^., = cosh0fcC^-sinh6'feC_fc 

where 

k k 

Unlike the Heisenberg ferromagnet, where the classical 
ground state is also the quantum mechanical one, quan- 
tum fluctuations are present for the XY ferromagnet. As 
a result (Pf) < 1/2 even for the ground state. Pro- 
vided the difference is small though, the classical ground 
state can be considered to be a good approximation and 
the spin wave theory should be reliable. As tempera- 
ture increases, however, the amount of fluctuations in- 
crease, and eventually (Pf ) calculated by spin wave the- 
ory becomes < 0, which clearly indicates spin wave the- 
ory ceases to be valid. However, spin wave theory can 
be applied if a Lagrangian multiplier or magnon chemi- 
cal potential A is added to restrict the total number of 
magnons excited. The modified spin wave Hamiltonian 
is given by 

H'sw = Hsw-XY.i^"''tci) (14) 
I 

with the provision that if {c^ cij < 1/2 when calculated 
without A, then A = 0, and if {c^ci) > 1/2, then A be- 
comes non-zero and is adjusted to enforce {cfci) — 1/2, 
i.e. {Pi) — 0. After Bogiiubov transformation, H'sw 
becomes diagonalized 

H's^ = -NZ^~\N + Y,EJdtdu+^^ (15) 



where the magnon excitation energy 




and A' = = I ^ e"^-^ . The Bogiiubov coefficients 

are found to be 

cosh20fe-^,sinh20fe = ^ (17) 

-C'fc i^k 

with 

Uk = hZTi^ - ^tbZfi'^jf. + A, = ^t^ZTi^^f. (18) 
The Bose condensed order parameter is 

(h) = {P') = l-^Y. co^h 20, {nk + ^ (19) 
where 

rik = {d+dk) = (exp (/3Sfe) - 1)"^ (20) 

is the Bose-Einstein distribution of the magnons, I3~^ = 
T (wc pick energy unit so that the Boltzmann constant 

fcs — 1). The dependence of (h^ on n as calculated 

from linearized spin wave theory with and without the 
Lagrange multiplier at a fixe temperature is shown in Fig. 

[5^. When the order parameter (l)^ is positive, denoting 

the SF phase, A = and the dispersion Ek represents 
a gapless mode and is linear in k at small wave number 
k. On the contrary, the excitation energy has a finite 

gap for nonzero A and hence for zero (h^ , indicating a 

MI phase. 

Before we set off to calculate the self-consistent pa- 
rameters B and n, a glance at the bosonic distribution 
rik alerts us that the spin wave method with chemical 
potential is still problematic. Since J (x n?, the effective 
spin temperature is T' — T/n^. Though the Lagrangian 
multiplier prescription is reliable up to T ^ J = 2t},fi^, 
it fails when T » J. However, the necessity to consider 
small n values can mean T » J even at low temper- 
ature. One major problem with the modified spin wave 
method discussed so far is that at high T', the nearest 
neighbor spin-spin correlation function B can actually 
take on the wrong sign. This problem has been recog- 
nized in literature, and a remedy has been proposed^^^^. 
The correct sign of B can be maintianed if when calculat- 
ing the spin correlation {P[P[^g ~^ Pi'Pi+s) ' '^^^ keeps 
not just the quadratic, but the quartic order in HP boson 



7 



0.4- (a) 




Normal Spin Wave 
Normal Spin Wave 
with Larangian Multiplier 
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Normal Spin Wave 




Normal Spin Wave 




with Larangian Multiplier 




Modified Spin Wave 
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FIG. 2: The modifications made in solving the spin XY model 
with 'vacancies', (a) The superfluid order parameter < b > 
and (b) the spin correlation B = {SfSj + SfS^) before and 
after modifcations as functions of n at kT — O.ObtbZ, where 
n is not yet self-consistent solutions. < b > is tuned to be 
non-negative at small n's by a Lagrangian multiplier. The 
correlation B is further corrected by the modified spin wave 
method. 



operators: 

{p{^p,%s + ptpr+s) (21) 



The quartic terms of c operators are then evaluated by 
Hartree-Fock (HF) decoupUng 

(22) 

In the theory by Takahasi, these HF terms need to be 
determined self-consistently, in essence renormalizing the 
quadratic Hamiltonian Hg^/. That is too complicated to 
do in our problem, so we just evaluate them using the 
unrenormalized H'gyy, in line with the scheme proposed 
by Hirsch and Tang2^. Using H'g^^, the HF terms can 
be derived as 



1. 
1 



(23) 



{c,Cj) = (c+c+) 



where 



k 



Tj)) cosh26'fc 



9{ri-rj) = y^exp {ik ■ [rj ~ rj)) sinh26'fc 



(24) 



Consequently, the spin correlation becomes 

B(n,T) (25) 

5-(/(o)44,(^)4/(^) 

+J(ff('5)-/(<5))' + ^(g(5)-^g(0) 
+^ (/W-^5(0))'-|.g(0)' 

Though B does not appear in the form of a complete 
square as the case of Heisenberg fcrromagnet'^^, the in- 
clusion of the quartic terms nevertheless ensures that it 
has the right sign at all effective spin temperature T' (See 
Fig. 12b)- Together with the equation for n Eq. this 
equation for B provide the self-equations for n and B. 
Finally, within our MFT, the free energy per site is 

G{T) = -—\nTrexp[-l3(Hn + C + Hs)] (26) 



N 



tbZn^B{n,T)-TSs 



-.) 



where Sg is the spin entropy. The free energy is use- 
ful in case where there is more than one MF solu- 
tions to determine which is the most stable solution. 
Within our MF approach this is derived in principle 
from H'gy^r, which gives another problem. At high tem- 
perature, the spin entropy for a spin 1/2 system satu- 
rates to In 2. However, using H'g^^ the spin entropy is 

Ss = X]fe{(l + "fe) 1^ (1 + "-fc) ~ "-fc lii'^fc} s-nd can exceed 
this saturated value due to large fluctuations even when 
the Legrange multiplier limits the average value of rxfe. 
Since the effective spin temperature T' can be enormous 
at small n, this problem is relevant when we need to 
compare the free energy of the n « to the n w 1 MF 
solutions at temperature of interest. Our remedy for this 
problem is to impose an upper limit Ss = A:Bln2 when 
the calculated value of S", exceeds that value. 



III. RESULTS AND DISCUSSIONS 

The results obtained for the 3D SC lattice are as fol- 
lows. The main result will be that i) the T = first or- 
der transition obtained by DMZ in their single-site MFT 
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remains robust with phase fluctuations; ii) close to the 
T = transition point, the finite temperature transi- 
tion will also be first order—, and iii) the stablization of 
the ordered phase gives rise to reentrance with increas- 
ing temperature into the SF phase on the MI side of the 
T = transition. We first present our results for bosons 
in optical lattice in which case the relevant free energy 
density to be minimized is the free energy per site with a 
fixed chemical potential /i. Within our model, the exper- 
imental change in the strength in the periodic potential 
corresponds to changing ^jty^Z. 

At zero temperature, there are no magnons. Using 
spin wave approximation, the superfluid order parame- 
ter and the spin correlation are respectively (&) ~ 0.478 
and B* = B {n,T = 0) ~ 0.265 for aU non-zero n's. The 
self-consistent Eq. (0) and its T = limit given by Eq. 
(HSl) indicate that when A/ttZ < B* , n ^ 1 will be the 
self-consistent solution that minimizes the total energy of 
the system, but when A/t^Z > B*, the ground state will 
have n — 0. However, the n = 1 self-consistent solution 
branch does not disappear and hence is metastable until 
A/tbZ increases past 2B* . Thus, there is a first order 
phase transition from the MI state to the SF state with 
decreasing A/ti,Z. Accompanying the MI-SF transition 
is a jump in boson density (or equivalently, vacancy den- 
sity) and exciton density . We thus recover the 'vacuum 
switching' first-order transition obtained by DMZ^ using 
the single-site MFT. Within single-site MFT, the tran- 
sition from MI to SF occurs at A/ti,Z = 0.25, which is 
less than our critical value of A/ti,Z = B*. This is the 
consequence of the SF phase being stabilized against the 
MI phase due to quantum phase fiuctuations. 

On the SF side, as the temperature rises from zero, 
thermal fiuctuations will decrease the SF order parameter 
(b) both directly and also indirectly through the decrease 
in n. At the same time, the decrease in spin correlation 
B can drive the system from the defect rich (n close to 
1) to the defect poor {n close to 0) phase. As T contin- 
ues to increase, we can expect one of two scenarios. The 
condensate amplitude may vanish while remaining in the 
defect rich phase, followed subsequently by a transition 
from the defect rich into the defect poor phase. Should 
that be the case, we expect (6) to go to continuously. 
The other possibility is that there is only one transition, 
which primarily is from defect rich to defect poor, and 
when that happens, it drives the Bose condensed ampli- 
tude to 0. Our calculation shows that it is the second 
case that is realized provided A/tbZ is not too far from 
B*. This is shown in Fig [31 The self-consistent solu- 
tions of n as T is increased are shown in Fig. [3Ja). We 
see that at finite T, in between the n values of the so- 
lutions that evolve from the n — and h — 1 solutions, 
there is a third self-consistent solution. This solution is 
unstable, corresponding to a local maximum in the free 
energy(Fig[3Kb)). The low and high n solutions are local 
minima in the free energy. At T = 0, the high n solution 
is the global minimum, but with increasing T, the small 
n solution becomes the free energy minimum at a critical 




-0.06- 

1 -0.08- 

-0.10 J, ■ , ' 1 ■- 1 ■ 1 ' r 

0.0 0.2 0.4 0.6 0.8 1.0 
<n> 

FIG. 3; Illustrations for A/hZ = S* - 2 x 10"^, where 
B* is the zero-temperature transition point. (a) The self- 
consistent criteria in which solutions are as the cross-over 
points of n (n) — cxp(-/3Bi) _ - temperature 

arises, a middle branch of self-consistent n appears. The 
low and high n branches correspond respectively to the nor- 
mal solid and supersolid solutions, (b) The free energy per 
site as a function of n. The normal solid and supersolid so- 
lutions have free energy around the two local minima. The 
first order SF-MI transition temperature is Tc/tbZ ~ 0.0246. 



temperature Tc. Beyond, Tc, the large n solution becomes 
metastable, but it remains a local minimum until a higher 
temperature where it merges with the with the unstable 
middle branch solution. The above behavior is typical of 
first order transition. In the present case, a jump in n 
occurs at Tc. Just below Tc, n is large enough that the 
coupling J in Hs is sufficient for (6) to be non-zero. Just 
above Tc. the jump to a small n results in a value of J 
such that Tc is above the Bose condensation temperature, 
and hence (b) jumps to across the transition. 

The value oi A/tbZ shown in Fig [3] is very close to -B*, 
yet Tc is not that close to 0. The reason for this is because 
of the presence of reentrance in the MI side of the transi- 
tion. The reentrance can be understood qualitatively as 
follows. Just on the MI side of the transition, the energy 
of the defect poor state and the defect rich state are al- 
most degenerate. The defect poor phase has excitations 
with gap A, and hence the difference of the free energy 
from the ground state energy is exponentially small at 
low T. The defect rich state however is a Bose conden- 
sate, and has gapless excitations. Thus, its free energy 
decreases by a power of the temperature. As the tem- 
perature is raised, the free energy gain can overcome the 
ground state energy difference between the two phases, 
and the defect rich phase becomes more stable than the 
defect poor phase. Thus as the temperature is raised, the 
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system first undergos a (first order) transition from the 
defect poor MI phase into the defect rich SF phase, and 
then later on at a higher temperature undergoes the SF- 
MI first order transition discussed previously back into 
the defect poor phase. 

On a more quantitative level, the free energy of the 
defect poor phase at low temperature is given by 

Go (T) = -fcsT (In (1 + 2 exp {-A/ksT))) 

On the other hand, for the defect rich phase, the domi- 
nant temperature correction is due to phase excitations 
about the condensate, giving 



Gi (T) w A - hZB* 



C(3) 



{kBTf 



where the Riemann zeta function C (3) = 1-202. If 
Gi(r) < Go{T), then the stable phase is the defect rich 
Bose condensed phase. On the SF side (A < tf^ZB*), 
Gi(0) < Go(0), and the difference is further enhanced 
at low temperature due to thermal fluctuations. On the 
MI side, Gi(0) > Go(0), and the MI is the stable ground 
state. But as T increases, the correction in Gi(T) 
dominate over the exponential correction in Go{T), and 
provided the ground state energy difference is not too 
big, Gi (T) becomes smaller than Go (T) beyond some low 
temperature, giving rise to the 'reentrance' phenomenon. 
As T continues to increase, the disordering effect of ther- 
mal fluctuations finally dominates so that there is a sec- 
ond transition back into the defect poor MI phase. The 
temperature dependence of the free energies of these two 
competing phases is illustrated in Fig. [4] The phase di- 
agram in the vicinity of the T — SF-MI transition is 
shown in Fig. [5l The reentrance on the MI side can be 
clearly seen. 

Although our results are based on our modified spin 
wave theory for the spin Hamiltonian H^, we believe the 
results are quite reliable. As a comparison, we consider 
the case oi ID, where Hs can be solved exactly using 
the well-known Jordan- Wigner transformation. The ID 
case by itself is also interesting because it can be realized 
experimentally in cold atoms. Using the Jordan- Wigner 
transformation, the spin 1/2 operators are mapped into 
spinless fermion operators 
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at 
9tm 



'St 



(27) 



where 



i-i 

i=l 



and = 5imA9h9m} = 0. Then the spin Hamil- 

tonian can be rewritten by the fermion g operators and 
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FIG. 4: Illustrations for A/ttZ = B* + 10 ^, where B* is the 
zero-temperature transition point, (a) The free energy per 
site plots of the normal solid and the supersolid self-consistent 
solutions as functions of temperature. There are two phase 
transitions, both first-ordered, as the temperature increases 
for this specific A. (b) The solved superfluid order parameter 
< fe > as a function of temperature, showing sudden jumps 
at the two first-ordered phase transition points. 



0.030 
0,025 - 




0.020 - 






0.015- 


Supersolid 




0,010- 




/ Normal Solid 


0.005 - 






0,000 - 








-3 -2 -1 


1 2 3 



Deviation of A / tZfrom B ,10 



FIG. 5: The finite-temperature phase diagram, obtained by 
comparing free energy per site, around the the zero tempera- 
ture transition point A/ti,Z = B* ~ 0.265. The tick labels of 
the horizontal axis denote the deviation in the unit of 10~^tbZ 
from the zero temperature transition point. There is a region 
on the zero-temperature normal solid side where the reen- 
trance phenomenon occurs at the finite temperature. 



be diagonalized in k space as 



Y^k9t9k 



with 



Ek — —tbZn cos k 



(28) 
(29) 



10 
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FIG. 6: The finite-temperature phase diagram around the 
the zero temperature transition point A/f;,^ = B* ~ 0.265, 
obtained by comparing free energy per particle rather than 
the value per site as shown in Fig. [5] The result compared 
to Fig. [5] shows moderate difference. 



The spin correlation is then found as 

^^^E'^o^^tofc) (30) 



where 



(fffe 5fc) = (cxp (/3efc) 



1) 



(31) 



is the fermion distribution function. At T = 0, the 
negative energy levels are filled, while the positive ones 
are empty . The zero temperature transition is at 
A/tbZ = B** = i « 0.318. We note that in addition 
to solving Hg exactly, this result is actually the exact re- 
sult for the Hamiltonian Eq. ([2]) . This is because [H, rii] 
= , so ni is a good quantum number for the ground 
state. Assuming no breaking of translational invariance, 
we then have either all ni = or all Ui = 1. Compar- 
ing the exact solution to the modified spin-wave method, 
which gives A/tbZ = B*jj = 0.314 as the transition point 
in ID, we see that the approximate modified spin wave 
method does very well. We also note that both meth- 
ods show that there is no Bose condensation in ID at 
T = 0, and the transition in n is driven by short range 
spin correlation B. 

At finite T, the free energy due to Hg can still be cal- 
culated exactly in ID, but the MFT procedure we use to 
decouple S and n is an approximation. The MF finite 
temperature results for the ID case by both the modi- 
fied spin wave and Jordan- Wigner transform are shown 
in Fig. [7] and [8l The two method gives similar pre- 
dictions in the spin correlation and the self-consistent n 
solutions. The results of free energy obtained by the two 
methods exhibit some difference, but the lower the tem- 
perature the smaller the discrepancy. More importantly, 
the reentrance behavior from defect poor to defect rich 
back to defect poor phases on the MI side appears in both 
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FIG. 7: ID MFT results of (a) correlation function B (n, T) , 
(b) the self-consistent scheme in finding n, (c) free energy 
per site and (d) free energy per particle by Jordan- Wigner 
transformation and the modified spin wave method. System 
parameters are A/tbZ = 0.3 and kT/tbZ = 0.1. 



methods. Both methods give phase transitions in n at fi- 
nite temperatures in ID, which is of course incorrect and 
an artifact of MFT, and the transitions should be inter- 
preted as crossovers in defect densities rather than real 
phase transitions. Nevertheless, the ID results provide 
additional support for the 3D phase diagram presented 
above. 

The results presented so far are relevant for bosons in 
optical lattice, where the experiments may be performed 
at constant chemical potential and the lattice constant is 
fixed by the optical lattice. To apply them to supersolid 
^He, modifications are necessary for the following reason. 
In the case of solid ''He, the lattice constant is not fixed 
externally but self-determined. Instead it is the number 
of helium atoms that is fixed. Therefore, when comparing 
the free energy of the defect poor and defect rich phase, 
one should use the free energy per atom rather than the 
free energy per site, with N — Na + Nb + Ny, where N is 
the number of sites, Na^b are the number of sites occupied 
by helium atoms in a and b states, and Ny the number of 
vacancies. The number of ^He atoms Nhc — Na + Nb- 
Minimizing the mean field free energy per atom instead 
of per site gives small quantitative corrections, but do 
not change the main results presented above. The phase 
diagram in that case is shown in Fig. [G] The transition 
temperature is changed, but the first order transitions 
and reentrance remain. 

In summary, we have shown that the MI-SF transition 
of a two-band boson Hubbard model can differ signifi- 
cantly from that of the one-band model. Instead of a 
continuous transition, the transition here is first order, 
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FIG. 8: ID MFT results of free energy per site of high n and 
low n self-consistent solutions indicating phase transitions by 
(a) modified spin wave approximation at A/tbZ — Bl^ + 
0.001 and by (b) Jordan- Wigner transformation at A/tbZ — 
B** +0.001. 



both for the T = and the finite T transitions. On the 
MI side, there can be a reentrant transition into the SF 
phase and then back into the MI phase as temperature is 
raised. The underlying physics is that the MI-SF tran- 
sition is actually a by-product of the transition between 
two competing phases, which are the defect poor and de- 
fect rich phases. In the defect poor phase, bosons are 
localized, while in the defect rich phase, they are delo- 
calized and Bose condense. The defect rich phase contain 
a finite density of excitons (defined in the Introduction 
section) and vacancies even at T = 0. Thus, the MI-SF 
transition is accompanied by a change in boson density 
in the case of optical lattice and a change from commen- 
surate to incommensurate solid in the case of solid ''He. 
Our results confirm that the results obtained by DMZ are 
robust when phase fluctuations are included, and also ex- 
tend their T = results to finite temperature. For 3D, 
these fluctuations arc taken into account using a modified 
spin wave theory. Although it is an approximate method, 
we believe it to be quantitatively reliable for our model 
by comparing the results in ID using this method and 
using Jordan- Wigner transformation. 

This work was partly supported by the Hong Kong's 
RGC grant 706206. 
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